% Function : Analyze Thermal Error data
% Paramter :
%   folder: folder path that contains data files(string)
% Return:
%   et: structure of Eddy and Temperature data(struct)
%   Fields of et:
%     Elabel: legends of eddy(cell array of string)
%     Tlabel: legends of time(cell array of string)
%     Time  : time data(vector)
%     Eddy  : eddy data(matrix)
%     Temp  : temperature data(matrix)
%     nEddy : number of eddy channels(double)
%     nTemp : number of temperature channels(double)
%     m     : number of samples(double)
%     pn    : 1 denotes that eddy is positively related to temperature,
%             while 0 is negative(logical)
function et = Analyze(folder, sw)
    id = 0;
    et = [];    % initialize struct et

    if ~exist(folder, 'dir')
        error(['folder ''', folder, ''' not exist']);
    end

    % check plot switch
    if nargin == 2 && strcmp(sw, 'plot')
        doplot = true;
    else
        doplot = false;
    end

    mpath = fullfile(folder, [folder, '.mat']);    % path of MAT file
    if exist(mpath, 'file')
        % show message
        id = id + 1;
        disp([num2str(id), ': Load struct et from MAT file']);

        load(mpath);
    end

    if isempty(et)
        % show message
        id = id + 1;
        disp([num2str(id), ': Contruct struct et from TXT file']);

        fpath    = cell(1, 4);
        fname    = cell(1, 4);
        fpath{1} = fullfile(folder, 'ET*.json');
        fpath{2} = fullfile(folder, 'ET*_Time.txt');
        fpath{3} = fullfile(folder, 'ET*_Eddy.txt');
        fpath{4} = fullfile(folder, 'ET*_Temp.txt');

        for i = 1 : 4
            fdir = dir(fpath{i});

            n = numel(fdir);
            if n == 0
                error([fpath{i}, ' do not exist']);
            elseif n > 1
                error([fpath{i}, ' more than one']);
            end

            fname{i} = fullfile(folder, fdir(1).name);
        end

        % read json file
        [et.Elabel, et.Tlabel] = getLegend(fname{1});

        et.nEddy = numel(et.Elabel);    % number of channels
        et.nTemp = numel(et.Tlabel);    % number of channels

        % read data files
        et.Time = load(fname{2});
        et.Eddy = load(fname{3}); et.Eddy = et.Eddy - et.Eddy(1, :);
        et.Temp = load(fname{4});

        [m1, ~ ] = size(et.Time);
        [m2, n2] = size(et.Eddy);
        [m3, n3] = size(et.Temp);

        % check whether length of samples is consistent
        if m1 ~= m2 || m1 ~= m3
            error('number of samples is not consistent');
        end

        et.m = m1;   % number of samples

        % check whether numbers of labels match data
        if n2 ~= et.nEddy
            error('number of Eddy channels do not match that of labels');
        end
        if n3 ~= et.nTemp
            error('number of Temp channels do not match that of labels');
        end
    end

    % plot eddy-time & temp-time graph
    if doplot
        % show message
        id = id + 1;
        disp([num2str(id), ': ETPlot']);

        f1 = ETPlot(et.Time, et.Eddy, et.Temp, et.Elabel, et.Tlabel);
    end

    % if there is no sect in et, then do section analysis
    if ~isfield(et, 'sect')
        % show message
        id = id + 1;
        disp([num2str(id), ': section analysis']);

        % section analysis
        % ======== deprecated ========
        % sect = cell(et.nTemp, 1); % raw sect
        % for i = 1 : et.nTemp
        %     sect{i} = Section(et.Temp(:, i), 0.2);
        % end
        % ============================
        l      = Section(et.Temp(:, 1), 0.2); % left endpoints
        r      = l;                           % right endpoints
        l(end) = [];
        r(1)   = [];

        et.sect   = [l', r'];
        et.nSect  = size(et.sect, 1);


        % plot each channel of Temperature and section lines
        % for i = 1 : et.nTemp
        %     figure(); hold on;

        %     plot(et.Temp(:, i));

        %     y  = ylim;
        %     ns = numel(sect{i});

        %     for j = 1 : ns
        %         x = [sect{i}(j), sect{i}(j)];
        %         plot(x, y, 'k--');
        %     end
        % end

        % merge sections
        % ======== deprecated ========
        % sect  = mergeS(sect);
        % nSect = size(sect, 1);
        % ============================

        % append monotonicity and gradient information to section
        % each row of section: index of left, index of right, monotonicity, gradient
        half = et.nTemp / 2;
        et.sect = repmat(et.sect, 1, 2);
        for i = 1 : et.nSect
            il   = et.sect(i, 1);   % index of left endpoint
            ir   = et.sect(i, 2);   % index of right endpoint
            grad = 60 * ((et.Temp(ir, :) - et.Temp(il, :)) / (et.Time(ir) - et.Time(il)));

            temp = sum(grad > 0);
            if temp > half
                et.sect(i, 3) = 1;          % monotone increasing
            elseif temp < half
                et.sect(i, 3) = -1;         % monotone decreasing
            else
                et.sect(i, 3) = 0;          % unknown monotonicity
            end

            et.sect(i, 4) = max(abs(grad)); % absolute gradient: degree per hour
        end

        % pn : whether each eddy channel is positive correlation with temperature or not
        % act: order of activity(how sensitive each temperature channel is)
        for i = 1 : et.nSect
            if et.sect(i, 3) > 0 && et.sect(i, 4) > 0.8
                il = et.sect(i, 1);
                ir = et.sect(i, 2);

                et.pn  = et.Eddy(ir, :) > et.Eddy(il, :);
                et.act = SortedIndex(abs(et.Temp(ir, :) - et.Temp(il, :)));
                break;
            end
        end
    end

    % search for spindle accurate-stop experiment sections
    % search for valid sections that are part of spindle accurate-stop experiment
    if ~isfield(et, 'expt')
        % show message
        id = id + 1;
        disp([num2str(id), ': search experiment sections']);

        lgc = false(1, et.nSect);
        for i = 1 : et.nSect
            if et.sect(i, 3) > 0 && et.sect(i, 4) > 0.35
                lgc(i) = true;
                if i < et.nSect
                    lgc(i + 1) = true;
                end
            end
        end

        if sum(lgc) == 0
            warning('no accurate-stop experiment');
            save(mpath, 'et');
            return;
        end

        i = 1;
        k = 0;
        et.expt = zeros(1, 2);  % spindle accurate-stop experiment sections
        while i <= et.nSect
            if lgc(i)
                j  = i + 1;
                while j <= et.nSect && lgc(j)
                    j = j + 1;
                end

                % find one section
                k = k + 1;
                et.expt(k, :) = [i, j - 1];

                i = j;  % i: fast move
            else
                i = i + 1;
            end
        end

        et.nExpt = size(et.expt, 1);  % number of experiments

        % delete useless data in the head and the tail
        for i = 1 : et.nExpt
            st = et.expt(i, 1);
            ed = et.expt(i, 2);
            % deal with start of section
            il = et.sect(st, 1);
            ir = et.sect(st, 2);
            st = et.sect(st, 1) + cutS(et.Temp(il:ir, et.act(1)), 0.8) - 1;
            % deal with end of section
            il = et.sect(ed, 1);
            ir = et.sect(ed, 2);
            ed = et.sect(ed, 2) - cutS(et.Temp(il:ir, et.act(1)), 0.6) + 1;

            et.expt(i, :) = [st, ed];
        end
    end

    % plot section results(base on ETPlot figure)
    if doplot
        figure(f1); hold on;

        subplot(2, 1, 1);
        y1 = ylim;  % y values of vertical line

        subplot(2, 1, 2);
        y2 = ylim;  % y values of vertical line

        for i = 1 : et.nExpt
            x = [et.expt(i, 1), et.expt(i, 1)];
            x = et.Time(x);
            subplot(2, 1, 1); plot(x, y1, 'b--');
            subplot(2, 1, 2); plot(x, y2, 'b--');

            x = [et.expt(i, 2), et.expt(i, 2)];
            x = et.Time(x);
            subplot(2, 1, 1); plot(x, y1, 'k--');
            subplot(2, 1, 2); plot(x, y2, 'k--');
        end
    end

    % if there is no gc in et, then do correlation analysis
    if ~isfield(et, 'gc')
        % show message
        id = id + 1;
        disp([num2str(id), ': correlation analysiss']);

        % correlation analysis
        et.gc  = zeros(et.nEddy, et.nTemp, et.nExpt);
        et.gco = zeros(et.nEddy, et.nTemp, et.nExpt);
        for i = 1 : et.nExpt
            % generate test data
            il   = et.expt(i, 1);
            ir   = et.expt(i, 2);
            Eddy = et.Eddy(il:ir, :);
            Temp = et.Temp(il:ir, :);
            for j = 1 : et.nEddy
                if ~et.pn(j)
                    Eddy(:, j) = -Eddy(:, j);
                end
            end

            % gc: gray correlation degree, gco: descending order of gc
            [et.gc(:, :, i), et.gco(:, :, i)] = GrayCorrelation(Eddy, Temp);

            % bar graph
            % figure; bar(et.gc(:, :, i));
            % xticklabels(et.Elabel);
            % legend(et.Tlabel);
        end
    end

    % save data to MAT file
    save(mpath, 'et');


    % regression analysis(write results to file)
    % show message
    id = id + 1;
    disp([num2str(id), ': regression analysis']);

    fid = fopen(fullfile(folder, [folder, '.txt']), 'w', 'n', 'UTF-8');

    for i = 1 : et.nExpt
        fprintf(fid, 'Experiment: %d\r\n', i);
        il = et.expt(i, 1);
        ir = et.expt(i, 2);
        X  = et.Temp(il:ir, :);
        for j = 1 : et.nEddy
            % print header
            fprintf(fid, '## %s ##\r\n', et.Elabel{j});

            Y = et.Eddy(il:ir, j);

            [I, R, RR] = CForward(Y, X, et.nTemp, et.gc(j, :));

            % print statistics
            for k = 1 : et.nTemp
                fprintf(fid, 'R=%f RR=%f P:%s\r\n', R(k), RR(k), et.Tlabel{I(k)});
            end

            % print endline
            fprintf(fid, '\r\n');
        end
    end

    fprintf(fid, '\r\n');
    fclose(fid);
end
